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Abstract 

The origin of the Finite Grid Instability (FGI) is studied by resolving the dynamics in the ID 
electrostatic Particle-In-Gell (PIG) model in the spectral domain at the single particle level and at 
the collective motion level. The spectral fidelity of the PIG model is contrasted with the underlying 
physical system or the gridless model. The systematic spectral phase and amplitude errors from 
the charge deposition and field interpolation are quantified for common particle shapes used in the 
PIG models. It is shown through such analysis and in simulations that the lack of spectral fidelity 
relative to the physical system due to the existence of aliased spatial modes is the major cause of 
the FGI in the PIG model. 
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1. Introduction 


N-body type problems arise in many disciplines and underpins our understanding of complex 
dynamical systems like plasmas and the cosmos. In a typical N-body problem, the interaction in 
particle pairs can be of electrostatic, or electromagnetic, or gravitational in nature and each particle 
responds to a force that is the linear superposition of all one-to-one interactions it receives. Direct 
calculation of all one-to-one interactions of Np particles has a computation cost of 0{Np), therefore 
it is amenable to numerical simulation only when Np is small. The PIC method [1], or more 
generally the particle-mesh method, is an efficient numerical method that reduces the computation 
complexity by introducing a computation grid and taking advantage of the linearity in the sum of 
one-to-one interactions. In the PIC method, the interaction among the particles is mediated by 
the grid through the Green’s function of the interaction represented on the grid. The computation 
complexity is reduced from 0{Np) to 0{Ng) + 0{Np)^ where Ng is the number of grid points. When 
the number of particles per cell Np/Ng 1, the gain in speedup is large Np), therefore the PIC 
method is a popular choice in the ab-initio numerical simulation of N-body systems. However, two 
major problems arise in the PIC method due to the discrete grid: (1) the use of an Eulerian grid 
for the moments of the particle distribution and fields, in conjunction with individual Lagrangian 
particles in continuous phase space, implies an inherent inconsistency; (2) the grid representation 
of the Green’s function is usually an approximation of the real Green’s function in the continuous 
space. 

Despite the computation efficiency of the PIC method and its wide-spread use, especially in 
plasma physics, common PIC models are vulnerable to an electrostatic numerical instability known 
as the Finite Grid Instability (FGI) [2l|3] (there is also an electromagnetic instability known as the 
numerical Cherenkov instability mEiE]). Early practitioners using the electrostatic PIC model to 
simulate plasma dynamics observed a heating effect to the plasma which depends on the numerical 
parameters, i.e., the grid size Ax and the number of particles per cell Nc- This numerical heating 
has been extensively studied since the early development of the PIC model and the empirical scaling 
of the heating time th of FGI in a thermal plasma, which has the form uOpTn rsj {\D/^xf{ND + Nc), 
is summarized in [1]. Here N^ is the number of particles in a Debye length oup is the plasma 
frequency. It is also known that FGI comes from the aliased modes in the system due to the 
incompatibility between the Fourier spectra of the discrete Eulerian and continuous Lagrangian 
variables. (Numerical Cherenkov instability may also come from an Eulerian-Lagrangian mismatch 
in the convective derivative [4]). The numerical instabilities have been conventionally analyzed 
as unphysical resonances between physical and aliased modes m n 0 El El E]. The locations 
and growth rates of the unstable modes have been solved for using linear dispersion analysis in 
limited, yet essential cases, i.e., for spatially uniform cold or Maxwellian distributions. It is worth 
noting that, unlike the two-stream type instabilities, for some commonly used electrostatic and 
electromagnetic PIC models, FGI can arise without the intersection of the physical and aliased 
modes BU- Analysis for more realistic and nonlinear simulations has not been carried out. 

Various methods and numerical schemes to mitigate FGI have been proposed, including intro¬ 
ducing grid interlacing m and random jiggling [inilll], employing higher order particle shapes 0 
and temporal/spatial filtering m, implicit time differencing and enforcing the energy conservation 
property of the numerical algorithm between time steps da EH ESI ES|. All these techniques have 
shown great promise with regard to the reduction of the instability growth rate. However, this is 
often achieved with substantial distortion/damping of the meaningful dynamics at the short wave¬ 
length scale or by sacrificing conservation properties (such as the loss of momentum conservation 
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in an energy conserving algorithm, which has long been debated in the development of the PIC 
models [1]). Recent rigorous work on energy conserving algorithm has led to a large improvement 
of the momentum conservation through nonlinear iterations m 

The above efforts notwithstanding, the important questions about how and where FGI arises 
exactly remain to be answered. Previous works treat aliased modes as inherent in the system and 
study their properties and corresponding mitigation method. However, the origin of the numerical 
instabilities is clearly unphysical. Therefore, in principle a simulation plasma should be contrasted 
with the underlying continuous system, which has the same number of particles and particle shape 
as the simulation plasma, to determine the origin of the numerical instabilities. We call latter the 
physical system in the following, as it obeys a Vlasov equation for the shaped particles, as long as 
the same particle shape is used in defining the charge density from the particle and the electric field 
on the particle. 

There are many choices about what to be contrasted between the PIC and the physical systems. 
Conservation properties and dispersion relationship have been used. It should be noted that one 
direct consequence of FGI in a PIC model without built-in energy (momentum) conservation prop¬ 
erty is, as can be expected, the gross violation of the energy (momentum) conservation. For this 
reason, recent efforts have been devoted to improve the energy and/or momentum conservation of 
the numerical scheme in order to control FGI. But it should be emphasized that conservation laws 
are desirable when eliminating the FGI, but they are neither necessary nor sufficient conditions. 
An isolated plasma system can exhibit various kinds of physical instabilities while strictly conserv¬ 
ing momentum and energy. Furthermore, as total energy is a global property of the underlying 
microscopic processes and only one constraint on the degrees of freedom (two if total momentum 
is also considered) in the phase space, to understand what gives rise to FGI and its consequences, 
we need a better resolution into the dynamics. Linear dispersion, in which the eigen modes with 
complex Fourier frequency can be viewed as a way to resolve the (linearized) dynamics, is a better 
choice. However, such analysis is limited to special cases of the particle distribution and small 
perturbation amplitude. Insight from such analysis for the improvement to the numerical scheme 
is useful but not easy to obtained and applied to more general situations. Recently, symplectic PIG 
codes [mnH] have also been developed, for which the symplectic structures of the Hamiltonian 
system is preserved. The symplectic structure may be a good choice to contrast the PIG and the 
physical systems, however, it is not clear how it is related to the numerical instabilities at present. 

In this paper, we will study the FGI in the ID electrostatic PIG models by spectrally resolving 
the dynamics at the single particle level, thus allowing us to identify the components in the model 
that lead to unphysical instability. The dynamics in PIG result from the superposition of the pair¬ 
wise interactions as in a physical system. The major components of an electrostatic PIG model — 
the charge deposition, the field interpolation and the particle pusher, all operate on a single particle, 
while the field solver can be viewed as operating on the spatial Fourier modes. Therefore the use of 
the particle and spectral resolutions are natural choices for this purpose. Such a representation of 
the PIG models is given in section and the spectral errors in PIG models are analyzed in section 
As an alternative to the individual particle representation, one can also choose the modes of the 
collective particle motion as a representation. This has the advantage that the plasma dynamics can 
then be viewed as the collective wave-particle interactions and such couplings in a physical system 
and in a PIG model can also be contrasted. We note that the deposition, field interpolation and 
field solver only involve spatial operations at a fixed time, while the particle update in the pusher 
is a temporal operation in continuous space whose stability and convergence can be verified to rule 
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out its role in FGI. To facilitate simulation comparison with the physical system, the gridless model 
miEniEi] is used, in which all components and elementary operations of the physical system are 
projected onto the finite Fourier basis. It is demonstrated in section that the lack of spectral 
fidelity in the deposition, field interpolation is the major cause of the FGI. Finally we summarize 
in section [5l 

2. Spectral representation of the PIC model 

2.1. Charge deposition 

Let’s first look at the charge deposition scheme in Fourier space to understand the effect of 
aliasing. We will see that the most important effect in Fourier space is the summation over all 
Brillouin zones which is the result of a convolution process between a continuous spectrum and 
a periodic spectrum over Brillouin zones. The sampling needed to go from continuous space to a 
discrete grid is the cause of the latter spectrum. 

In a grid-based model like PIC, the contribution of a particle at position iTp^and of total charge 
Q on the density grid Vp is p{rp) = QW{rp^ Xp)^ where Vp is a vector on a uniform grid with grid 
size Ax = (Ax, A^, Az). The interpolation function for the deposition is W{rp^ r) — W{\rp — r\). 
Note that we have assumed the cell volume V — AxAyAz — 1 and dropped it for clarity. The 
difference between a particle shape function S{r) and a interpolation function W{rp^ r) is discussed 
, In the rest of this section, we will not distinguish these two and will use S{r) for 
clarity, p{rp) — QSp{rp — Xp). We define a transform (note this is not necessarily the proper 
Discrete Fourier Transform as Vp may be shifted from the origin of the coordinate system, as will 
be shown later). 


m 


Appendix A 


m = = QY,Sp{rp - x. 


)Q-ik-r 


which can be viewed as the continuous Fourier transform of p{r)^5{r — rp)^ where p{r) 
Xp). Then the non-unitary inverse transform is 


( 2 . 1 ) 

QSp{r- 


1 


'pdk. 


As Sp{r) is a function in continuous space. 


1 


( 2 . 2 ) 


Then, 


p{k) = ^ 5 ]( • dk'. 

oo oo 

(2.3) 


Carrying out the infinite sum (see Appendix B), we obtain 


^Italic bold font is used for vectors. 
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( 2 . 4 ) 


m = Qj2i^ici-kg, Arp)S,{kg)e-^’^‘>-^^, 


and 

where 


p{k + q' • kg) = V’(-q' • kg, Arp)p{k), 
V’(q-feg, Arp) = 


(2.5) 


is a phase factor. Vp = Vf + Avp and is a reference grid that contains the origin of coordinate 
system, kq = fc + q-fc^ is an vector alias of k in the q-th Brillouin zone [q-fc^ —fc^/2, q-fc^ + fc^/2), 
kg — (27r/Ax, 27tIA y^ 27vlAz) and 


n 


q = 


m 


I 


In the case that Avp ^ 0, the spectra in Brillouin zones are not necessarily the same as a result of 


the transform defined in Eq. (2.1). 


For a physical system consisting of particles of shape *Sp(r), a particle at position Xp will 
contribute to the density as, 


p{k) = QSp{k)e 


— ikXr: 


(2.6) 


where k G (—oo, oo). For a gridless model which does charge deposition in Fourier space using 


truncated Fourier basis, Eq. (2.6) is used for a truncated domain k G [— fcc? where kc is the 


cut-off wavenumber. To compare with a grid-based model, it is fair to set kc = kg/2 so both 
models have the same spectral resolution. 


In Eq. (2.4) for the grid-based model and Eq. (2.6) for the physical or the gridless model, one 
can clearly see that the charge density is from the contribution of a point-charge particle 

modified by its shape S{k). Furthermore, the grid-based model has the aliasing effect in Fourier 
space which is absent in the physical or the gridless model, i.e., the summation over all Brillouin 
zones. The Brillouin zones exist due to the discreteness of the grid. The summation is due to the 
need to convolve the particle shape’s continuous spectrum with a spectrum that includes Brillouin 
zones. 


The ratio between p{k) for the grid-based model, Eq. (2.4), and gridless models, Eq. (2.6), is 


(5p(fc))-ie*^--^V’(q- fcg, Arp)Sp{kg)e-^^^-^^ = 

q q 

= G{xp - Arp, fe; Spik)). (2.7) 

G{x^ k] Sp{k)) determines the spectral error in the deposition for a particular particle shape Sp{k) 
and will be quantified in section 

Eqs. (2.4) and ( |2.6[ ) give the density for a single particle. For an infinite number of particles 
from a distribution /(xp,i;), we have. 


5 









p{k) = f dvdXpf{Xp,V)Qj2^{q-kg, Arp)Sp{kg)e-^'^0-^^ = f dvQj^i^id-kg, Arg)Sp{kg)f{kg,V) 

q q 

( 2 . 8 ) 


for the PIC model and, 


p{k) = j 


dvQSp{k)f{k,v) 


(2.9) 


for the physical or gridless model. Apparently Eq. ( |2.8[ ) and (2.9) differ in both amplitude and 
phase, however, such differences depend on f{k^v) and are difficult to quantify. 

2.2. Field solver 

For a physical system, Poisson’s equation and electric field in continuous space are 

Vr : -VV(r) = p(r), E{r) = -V^{r). 

In Fourier space, this is completely equivalent to. 


Vfc : |[fc]r0(fc) =p(fc), E{k) = -i[fc]0(fc). 


( 2 . 10 ) 


where the operator [fc] = fc. 

In a simulation, there are many ways to solve for the field. The common choices are, 

• Use a spectral solver, 

Vfc G [-kg/2, kg/2] : \[k]sp\^{k) = p{k), E{k) = -i[k]sp^{k) (2.11) 

where [k]sp = k. In this case, E and p are collocated on the same grid which is chosen to be the 
reference grid, i.e., ve — Tp — ry, E{k) — '^E{rf)e~'^^''^f and p{k) — ^p{rf)e~^^''^f. We also 

have E{k + kg) — E{k) and p{k + kg) — p{k). 

• Use finite difference on Poisson’s equation and the potential for which the spectral represen¬ 
tation is, 

Vfc : \[k]fd\^{k) = p{k), E{k) = -i[fe]/d^(fc). (2.12) 


Here, [k\fd = {\kx\fd, [ky\fd, [kz\fd), [ka]fd for a G {x, y, z} is related to the Fourier representa¬ 
tion of the finite difference operator. p{k) — '^p{rp)e~'^^''^p and Ea{k) = 


PEa 


from the transform defined in Eq. (2.1), which in turn gives the usual form of [ka]fd^ o.g., 
[kx]fd — kxSmc{kxAx/2). As central differencing is typically used, [kx]fd^ [^y]fd: [f^z]fd ^^re purely 
real functions. Since Ea and p may be defined on different grids from the reference grid, so in 


general Ec{kFkg) = ^{-kg, ArEjEa{k) , p{k + kg) = ^{-kg, Arp)p{k) and Ave^ = ^E^-Vf, 

Avp = rp- Vf. 

Another possible choice is to use spectral solver on Poisson’s equation and finite difference on 
the potential, as done in the code XESl [1]. It is clear that none of these choices for the field 
solver has systematic phase error, however, using finite difference inevitably introduces systematic 
amplitude error. 
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2.3. Field interpolation 

In PIC, after the field on the grid E{rE) is calculated, it is interpolated to the j-th particle’s 
position according to 


= F.SEixj - rE)E{rE), (2.13) 

ve 

where the particle shape function is S'^(r) = f^^dkSE{k)e'^^''^/S tt^ and may be different than 
Sp{r) in charge deposition. Then, 

I roo -I POO 

= dkSE{k)e^^'^^^E{rE)e-^’^'^^ = ^ j dkSE{k)e^’^-^i E{k). (2.14) 

This means that S{k) — SE{k)E{k), i.e., no aliasing effect for a particular mode E{k). This 
is consistent with the idea that only sampling leads to aliasing, interpolation does not. Note that 
this view is for single mode analysis, for the total force on a particle, we need to integrate over 
all fc, so the summation over all Brillouin zones reappears. Generally, we have Ea{k + kg) = 
^{—kg^ ArEc)Ea{k) for each component a = z, so 

Saixj) = ^ dfe^«(fc) VV’(-q • kg, ArEjSEikg)e^^’>-^i. (2.15) 

8vr J-kg/2 “ 

In a physical model, the field on the particle can be directly reconstructed from the Fourier 
components. 


1 

^(*i) = ^ y dkSE{k)e^'"'^^E{k). 


(2.16) 


In a gridless model, only modes from the truncated Fourier basis, i.e., k G [—kgl2^ are 

used. 


^(®i) = A dkSE{k)e^'"-^iE{k). (2.17) 

J-kg/2 

Similar to the charge deposition, the integrands in the above equations for the PIC, physical 
(or gridless) model differ by a ratio, 

{SE{k))-^e-^^-^iY.'^{-ci-kg, ArEySE{kg)e^^’>-^i = 

q q 

= G(Ar£„ - Xj, k; SE{k)). (2.18) 


3. Spatial phase and amplitude errors in PIC 

Depending on the initial conditions, driving force, equilibrium etc., the plasma that we are 
trying to model may exhibit a host of unstable modes itself. So to distinguish the physical and 
numerical instabilities, one needs to look for the differences between the two systems. As the spatial 
components of the PIC loop, i.e., charge deposition, field solver and field interpolation, are now 
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spectrally resolved, it is informative to investigate the errors in the spectral domain involved in 
these components. In fact, it turns out that all these three components have errors in the spectral 
domain. Such errors can be systematic or random in nature (e.g. from round-off error) and we will 
only focus on the former. 

It is particularly worthwhile to understand the role that the systematic spectral errors play in 
the numerical instabilities. The phase error plays a crucial role in the development of instabilities, 
and in many cases determines whether the system is stable or not. On the other hand, the amplitude 
error typically affects the instability growth rate. In some other cases, e.g., for instabilities with 
amplitude threshold or some parametric instabilities, the mode amplitude can also determine the 
stability. However, most numerical instabilities we encounter and need to mitigate in the PIC model 
are not known to be of this type. 


3.1. Spectral error for shaped particle 

Next, we will quantify the phase and amplitude errors of the PIC model compared to the gridless 
model, which is a computationally feasible but costly model for approximating the physical system 
with both spectral fidelity (up to the cut-off wave-number) and conservative properties. 

For the field solver, it is easy to see from Eqs. ( 2.11[ ), ( |2.12 ) that neither the finite difference or 
spectral solver introduces systematic phase error. There is systematic amplitude error in the finite 
difference solver but none in the spectral solver. 

In general, Eq. (2.8) indicates that the deposited charge density has systematic amplitude 
and phase errors which depend on the particle distribution function. But for a system with finite 
number of shaped particles, we can further discern these errors by resolving the contributions from 
individual particles using Eq. (2.7). As summing each individual particle’s contribution to the 


charge density is a linear process, we do not expect such summation to generate additional errors 
other than round-off. In addition, field interpolation is naturally written for one particle in Eq. 


(2.15). Therefore, for a single particle, Eqs. (2.7) and (2.18) give the errors that can be quantified for 


a specific particle shape. For an orthogonal coordinate system and corresponding deposition and 
force interpolation, G(x, fc; S{k)) = Yl Sa{ka)) for S{k) = Yl For 

a=x,y,z a=x,y,z 

Cartesian coordinate, we only need to determine the spectral errors for one particular dimension, 

e.g., g{xlAx^ k] s{k)) = for x direction. 

q 

Since g{xlAx^ k] s{k)) is periodic in x, 

g{x -h 1, k] s{k)) = g{x, k] s{k)), 
and symmetric in k for symmetric shape function s{—k) — s{k), 


g{x, -k; s{-k)) = g{x, k; s{k)), 

one only needs to determine g{x^ k] s{k)) for x G [0,1] and k G [O^kg/2]. For the most commonly 
used B-spline particle shapes Sg{k) of order g and the Gaussian function sg{k)^ the systematic 
spectral errors g{x^ fc; s{k)) can be analytically quantified. 

For B-spline particle shapes, Sg{k) = [sinc(A:Ax/2)]^+^, the infinite sum in ^(x, k; s^) can be 
evaluated into a compact form using the Hurwitz-Lerch transcendent $( 2 ;, s, a) = z'^{n + a)~^ 

to give. 













g{x, k, 5^) = ^+ 1,1 + 


+$((-ir+ie*2-,^ + l, 


(3.1) 


Similarly, for a Gaussian function of width a Ax, sg{x;a) = sg{k;a) = 

^-{kaAxl2) the infinite sum in g{x, fc; 5 ^) can be evaluated using the elliptic theta function 

g) = 1 + 2 cos{2nu) to give, 


g{x, k, sg) = (kAx/2 + ix/a‘^, . 


(3.2) 


Although sg{x] a) is a valid particle shape but not a valid interpolation function, one can define a 
new function s{k) — sg{k] (j)s^{k) such that the corresponding s{x) is a valid interpolation function. 
For example, for /i = 0, we have SErf{k]^) = sg{k](j)s^{k) — e-(^cr^^/ 2 ) sinc(A:Ax/2)/\/^ and 
SErf{x-,cr) = [Erf(^ - j^) + Erf(^ + 4 ^)] /{2Ax) , where Erf(a:) is the error function. It is 
not clear whether g{x^ k; SErf) has a compact analytic form but numerical evaluation in Fig. 
shows g{x, k] SErf) is similar to g{x, fc; 5^). 

The amplitude and phase of g{x, k, 5^^), g{x, k, sg) and g{x, k] SErf) shown in Figs, [^and 

The general trends in Figs. and are that (1) larger amplitude and phase errors occur at 
higher k] and (2) the wider the particle shape, the errors are more concentrated at high k. Both 
amplitude and phase errors depend on the particle’s position in the cell, except for the amplitude 
error for sq (nearest grid point). For the phase error, the largest error occurs when the particle is 
in the middle of the cell. 

Comparing g{x, k, Sjj,) in Fig. another salient observation is that the phase error improves 
slowly with the increase of the order of the particle shape. Also the dependence of the errors on 
the particle position does not seem to go away quickly as the order is increased, despite stronger 
attenuation to the high k modes. As will be shown at the end of section |3.2| the spectral errors, 
especially their dependences on the particle position, play a key role in determining the growth rate 
of most unstable mode in the linear FGI ca se, w hile t he at tenuation to the high k modes due to 
S{k) and S'{k) in the front factors in Eq. (2.7) and (2.18) enlarges the stability domain, but at 
the expense of accuracy of the high k modes. The high order B-spline particle shapes offer both 
benefits but the latter benefit can also be achieved via a low-pass filter that is independent of the 
particle position. 

3.2. Relation between the linear instability and the speetral error 

The value of s{k) in and outside the fundamental Brillouin zone —kgl2 < k < kgl2 play different 
roles in the linear grid instability. The latter value is related to the spectral errors g{x, k; s{k)). In 


order to distinguish their roles, it is instructive to inspect the linear dispersion equation Eq. (C.7) 
in the limit At ^ 0 for a ID cold drifting plasma with velocity vq, 




'(/Cg)fc2 


- kqVoY 


= 0 , 


(3.3) 


where we also assume Sp = se = s, [kf — k^ and K{kq) = kq for simplicity. 
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|g(x, k, so)| 



n 1 
|g(x, k, si)\ 



n 1 

|g (X, k, S2)l 



n 1 
|g (x, k, S3)| 






Figure 1: The amplitude (left column) and phase (right column) of g{x, k, s^) for /x = 0, 1, 2, 3, which correspond 
to nearest grid point, linear, quadratic, cubic particle shapes, respectively. 
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Figure 2: The amplitude (left column) and phase (right column) of g{x, k; sg) for the Gaussian shape sg{x;Ax) 
(top row), sg{x;2Ax) (middle row) and of g{x, k; SErf) for the SErf{x;2Ax) shape (bottom row), respectively. 
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As = s{k)g{x/Ax,k]s{k))^ inverse Fourier transform with respect to x gives 

q 

s{kq) = gqs{k)^ where ga = k] s{k)) = Ax f dxe~'^^^^^g(x, k] s{k)) is discrete at cr = g due to 
the periodicity of g{x^ k; s{k)) in x. For nth order B-spline shape, gq — {—k/kq)'^^^] for S Erfj k] cr) 
shape, gq = —k/kq • e^^[{k^x/2)^-{kAx/2^q7r)^] ^ Keeping only the 0th and qth terms, Eq. (3.3) can 
be rewritten as. 


1 


\k) 


Zlkl/e 


1 

^ - kgvof 


= 0 , 


(3.4) 


where Q. — uj — kv{) and —kgl2 < k < kg/2. 

In Eq. ( |3.4| ), s{k) and gq affect the instability domain and growth rate differently. For the model 
dispersion in Eq. (3.4), mode stability requires s{k) < kgVo (see Appendix C). So for sufficiently 


small s{k) in the fundamental Brillouin zone, there is no instability for a particular k. However, 
this is achieved at the expense of the dispersion accuracy of the physical modes, since s{k) can 
be seen as a renormalization of Q. While gq = s{kq)/s{k), which is closely related to the spectral 
errors g{x^ fc; s{k))^ affects the alias modes thus the instability only. More specifically, gq has a 
large effect on the instability growth rate but little on the stability domain, gq is real and can be 
written as a sum of the Fourier coefficients of the amplitude and phase components of g{x^ k; 5(fc)), 


Sq = '^■^l^q-l 
I 


(3.5) 


where A/ = Ax f dxe fc; 5(A:))| and &q-i = Ax f dxe l)x-\-iATg[g{x,k; s{k))] ^ 

dominating terms in Eq. (3.5) is AqO^ and Ag©o as both A/ and ©/ are highly peaked at / = 0. 


Therefore the qth alias mode appears mainly because the presence of the amplitude or phase errors 
that is position-dependent and modulated at the qth harmonics of the grid. For the common B- 
spline particle shapes, Ao©^ dominates at high /c, therefore the amplitude error is more important, 
while at lower A:, A^©o dominates, i.e., the phase error is more important. 

The most unstable mode k^^^ in the dispersion Eq. (3.4) can be obtained from s{k'^^^) kgVQ. 
Its normalized growth rate is 


31/2 


, .max ^ 






,2/3 

>9 


31/2 


{kgvk)^'^ (1 + qkg/k^--f^ g2/3. (3.6) 


Therefore, smaller s{k) in the fundamental Brillouin zone (but keeping gq the same) actually leads 
to smaller k'^^^ and higher maximum growth rate, while smaller |g^| (but keeping s{k) in the 
fundamental Brillouin zone the same) lowers the maximum growth rate. This indicates that a 
particle shape with a high contrast between the fundamental and higher order Brillouin zones, i.e., 
small |gg| that translates to small position dependent spectral errors, is beneficial when aliasing is 
unavoidable. 


3.3. Spectral fidelity of the collective modes 

Now it is clear that the systematic spectral errors in the deposition and field interpolation 
occur at the individual particle level, but the overall errors from an ensemble of plasma particles 
are less transparent. To this end, we can analyze how these errors manifest themselves in the 
collective motion of a finite number of particles. To simplify this analysis, an initially cold beam 
of Np particles moving at velocity Vq is assumed, but the analysis is not restricted to this velocity 
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distribution. Each particle has charge Q, mass m. The system is periodic with length L and Ng 
cells, Ax = L/Ng is the cell size. When there is no perturbation, the j-th particle is located at Xj 
and particles are equally spaced with distance AX = L/Np. The particle’s position xj is perturbed 
at t = 0 and then evolves self-consistently. The displacement of the particle is 6x{Xj^t) = Xj{t) —Xj 
and can be expanded into Fourier modes 

5x{Xj,t) = AT-i S'x{K, = 2N-^ Y ^)) > 

where Sx{K,t) = = ^Sx{Xj,t)e~^^^A and 6{K,t) are the mode amplitude 

3 

and phase of the collective motion of these Np particles. The drift motion of the beam is described 
by the X = 0 mode. 

Assuming these particles have shape s{k)^ their contribution to the density of mode k in the 
periodic physical system (or a gridless model) is 


p{k) = qYk^y 


—ikxj{t) _ 


qk^)Y^ 


—ik 


Xj+2N-^ 


E A(K,t)cosiKXj+eiK,t)) 


(3.7) 


where |A:| < tv I Ax. the phase factor in this sum can be written, using the Jacobi-Anger expansion, 
as 


—ik 


Xj+2N-^ X; A(K,t)cos(KXj+0(K,t)) 


_ ^-ikXj 




. AvK[KXj+0(K,t)] 


.i)]| 


V0,...,VK 


JY:{^KK)-Xj+iY,UK-[0{K,t)+7T/2\ . 


n J,^{-2kA{K,t)/Np 


Uo,...,UK 


(3.8) 


where vk is an integer for each mode K in ^ X ^ 0, Jy^{x) is the z/^-th Bessel function of 
the first kind. ( H ) is a sum (product) over all combination of z/^s and ^ z/^ is a sum 


1A0,...,1AK 

of every z/^ foi* 7^ > X > 0. 


Uo,...,UK 
TT 

AX. 


Summing Eq. (3.8) over Xj and using ^ 


^'0=—00 


Xo (a;) = 1 gives, 


P{k) = NpQsik) Y J](i/kK)) . Y[ J,^{-2kA{K,t)/Np) 

(3.9) 

where kp = 2tv jAX — Np2TvjL and m is an integer. 

Eq. ( |3.9[ ) describes the contribution from the Fourier modes of the collective motion to the 
Fourier modes of the density, where the first, second and third terms in the curly bracket determining 
the mode coupling relation, phase and strength, respectively. The 7r/2 in the phase factor comes 
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from the phase difference between particle’s displacement and the density. Note that Eq. (3.9) is 


valid for arbitrary number of particles and arbitrary spectrum/amplitude of the collective modes. 
Due to the finite number of particles in a system, the mkp term in the argument of the delta 
function introduces a physical coupling due to the “particle aliasing” effect of the collective mode. 
Typically we are interested in the case /c <C /cp, thus this effect couples a low k mode and a very 
high K collective mode efficiently only when z/p is a small number. Furthermore, the coupling 
coefficient {—2kA{K^t)/Np) would be small when there are sufficient large number of particles, 
Np ^ so practically we may take m — 0 and ignore this effect in a physical system for most 
situations. 

Next we investigate the coupling among collective modes and density modes only. This can 
be illustrated by understanding how a single collective mode couples to the density perturbation 
first. For a collective mode spectrum with a single mode Kq^ A{Ko) ^ 0, while A{K ^ Kq) = 0. 
Since ^ 0 only when z/p = 0, so we have S {k — ^ (ukK)) — 5 [k — This cor¬ 

responds to the excitation of the fundamental mode k = Kq and its harmonics k = z/p^i^o 
with the mode coupling coefficient being s(z/p:o^o)e^^^o * J^kq {~‘^^KoKoA{KoA)/^p)- 

Similarly, when there are two collective modes Kq and i^i, their beating will result in den¬ 
sity perturbation at mode k — pkqKq + i^KiKi with the coupling coefficient being s{v'KoKo + 

{-2kA{Ko,t)/Np)J^j^^ {-2kA{Ki,t)/Np). From 
this, the coupling for the density perturbation involving more collective modes can be generalized. 

The density perturbation created by the collective modes in turn drives more collective modes 
of the particles through the electric field. We define the electric field on the particle as a function 
of X^-, 


whose Fourier transform with regard to Xj is. 


S{K,t) = Ng 


J[kxj{t)-KXj] 


Here, [k] is the effective operator from the field solver. Similar to the phase factor for the density 
in Eq. ( |3.8[ ), the phase factor in the above equation can be expanded and the sum over Xj can be 
carried out to give. 


= E \s{K-k-mkp-J2{uK'K' 


'[k] 






n 


{2kA{K',t)/Np) 


(3.10) 




The relation between S{K^ t) and p{k) (or E{k)) depends on the argument of the delta function 
and the sum over A;, therefore it is more complicated than that between p{k) and 6x{KA) in 
Eq. ( |3.9[ ). Again we may ignore the particle aliasing effect in the collective mode when N ^ 1. 
Then for the simplest case of a collective motion with a single mode Xq, a non-zero value of 
n {‘2kA{K'A)/Np) requires that z/^y^o = 0 ^ therefore, J2{^k'K') = 
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^KqKq. Coupling to £{K^t) thus comes from p{k = K — i^KoKo)^ which is not necessarily the 
harmonics of Kq. When the collective motion involves two modes Kq and coupling to £{K^t) 
can come from p{k = K — i^KoKo — ^KiKi) for any vkq and vki- 

Finally, the evolution of the collective modes is determined by the equation of motion 


didxiK^t)) / dt^ — ^£(K^t)^ 
m 


(3.11) 


which is an ordinary differential equation (ODE). Clearly, regardless of whether it is written in a 
continuous or discrete time variable, this step does not introduce spectral errors in iF-space. 

Now, Eqs. (3.9), (3.10) and (3.11) are the coupled equations describing the evolution of the 


collective modes in the physical system of a cold beam. When the discrete time version of Eq. (3.11) 
is used, these equations represent a discrete-time analog of the physical system. If Eq. (3.11) is 


discretized using a scheme consistent with the underlying ODE, i.e., the solution of the discretized 
equation approaches that of the ODE when the discretization step At ^ 0, then a numerical 
instability, if occurring in this discrete-time system, will exhibit dependence on At. Therefore, a 
convergence test can discern a numerical instability due to discretization in time. Early study of 
FGI already established that the numerical heating rate is independent of the time step, therefore 
time discretization cannot be the reason of FGI. From the single particle analysis of the systematic 
spectral errors presented earlier in this section, we can see that another possibility of the numerical 
instability in PIG is the spectral errors in the deposition and field interpolation due to the aliasing 


effect. In fact, for the PIG model, Eqs. (3.9) and (3.10) take on the following forms. 


p{k) = NpQ^s{kq) |(5 (fcg - mfcp - 

■ n J,^{-2kqA{K,t)lN^ 


[e{K,t)+TT/2] 


(3.12) 


Ul,...,UK 


N„ 


= /E E T^fikg)-s{kq)- \s(K-kq-mkp-J2{^K'K' 

^ q \k\<kg/2^ 






n J,^,{2kqA{K',t)/Np)\, 


(3.13) 




where k has been replaced by its alias kq and a sum of all aliases is taken. 


Gompared with Eqs. (3.9) and (3.10), Eqs. (3.12) and (3.13) lack the spectral fidelity to the 


physical system, i.e., they contain not only phase and amplitude errors, but also unphysical mode 
couplings. As will be shown in the next section, such lack of spectral fidelity is the major cause of 
the finite grid instability. 


4. Comparison of the gridless model and the PIC models 

In this section, we compare the ID non-relativistic electrostatic simulation results of a numerical 
plasma that is susceptible to the finite grid instability, for various PIC models and from the gridless 
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model. The simplest simulation is that of a single-mode electrostatic oscillation in a cold plasma 
with a drift. While this mode will generate high harmonics, such system should be stable for a long 
time. During this period, the plasma should remain cold, i.e., the local momentum spread is zero. 
The three spatial components of the numerical scheme — the charge deposition, the field solver 
and the field interpolation, are implemented in the PIC/finite-difference mode and/or the gridless 
mode. The mode history from the simulations are compared, the focus here is not the accuracy 
of the simulation but to understand the roles of the amplitude and phase errors in the finite grid 
instability. 

The simulation box has dimensionless size L — and N = 33 cells. The initial plasma 

density no is constant across the box. Nq = 300 particles per cell is used for the electrons and 
the ions constitute a uniform neutralizing background. The electrons are equally spaced and a 
sinusoidal displacement 6x{xpo) = ALcos{27rMxpo/L)/{27vM) is applied to initialize a density 
perturbation, where M is the mode number and is the unperturbed particle position. The 
density perturbation is given by ni(x)/no = —ddxjdx — Asi]i{2'kMx/L) for point particles and 
the actual amplitude will be smaller by s{k — 2 'kM/L) when a shaped particle is used in the 
simulation. The electron initial velocity is the sum of the drift velocity Vq = 0.01c and perturbed 
velocity 5V = ALujpSm{2'KMx/L)/{2 'kM)^ where c is the speed of light. 

Results are compared in Fig. for four simulations with UpAt = 0.2: (1) the momentum 
conserving PIC model [1] with linear particle shape; (2) the momentum conserving PIC model 
with quadratic particle shape; (3) the energy conserving PIC model [13] with quadratic particle 
shape; and (4) the electrostatic gridless model [20l |2T] with quadratic particle shape. For all PIC 
models compared, the initial small-amplitude single-mode perturbation {A 10 ^ for mode K — 9) 
excites all the modes in the system at a level of ^ 10“^ immediately after the first step. Then, 
it is observed that mode /c = 15, which is an alias mode from the second harmonics K — 18 of 
the initial perturbation in the adjacent Brillouin zone, grows exponentially in all PIC models. The 
growth rate is higher for the lower order particle shape used in the momentum conserving PIC 
models. Although the initial amplitude of mode /c = 15 is higher in the energy conserving PIC 
model with quadratic particle shape due to the lower order, i.e., linear, particle shape used in 
field interpolation, the growth rate is about the same as in the momentum conserving PIC with 
quadratic shape. Following the growth of mode k = 15, mode k — Q also grows at a similar rate, 
which appears to be a result of the beating between modes A: = 9 and k — lb. Other exponential 
growing modes before uopt < 50, e.g., modes k = 3 and A: = 12 may also result from beating. These 
modes grow to an amplitude comparable to the initial perturbation during ujpt ^ 50 — 100, after 
which a wide spectrum of modes are excited. As a result of the interaction with a wide spectrum 
of waves, the phase space plots in Fig. show heated electron distributions in all PIC models, 
while the momentum (energy) is relatively constant in the momentum (energy) conserving model 
by design. Comparing with the gridless model for which the perturbation remains single-mode and 
the electrons stay cold during the simulation, it can be seen that the development of the wide-band 
spectrum and the heating of the electrons are clearly numerical effects. Such numerical artifacts 
result from the lack of spectral fidelity in the PIC model. 

It should be noted that the gridless model improves not only the spectral accuracy of the 
dynamics, but also the energy and momentum conservation of the simulation model. In Fig. the 
gross loss of energy (momentum) conservation in the momentum (energy) conserving PIC models 
are correlated with the saturation of the unstable mode growth. Furthermore, despite the 2^^ 
order accurate leap-frog scheme used for the time advance in all models, only the gridless model 
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exhibits energy conservation (and possibly for momentum conservation) that scales like O(At^) 
in this numerical instability test. For the gridless model, the energy conservation is improved by 
about a factor of 1000 when the time step is decreased by a factor of 10 to cjpAt = 0.02, while 
the momentum conservation is already at the machine precision for ojpAt = 0.2 thus does not 
improve further. For the energy conserving model, the heating in phase space, mode history and 
growth rate are essentially unchanged when the time step is decreased to ojpAt = 0.02, while the 
momentum conservation is only slightly improved and the energy conservation is improved by about 
a factor of 10. The heating and unstable mode growth behave similarly in the momentum conserving 
model, but neither the energy or momentum conservation improves as time step is decreased to 
(jjpAt = 0.02. The superior conservative properties of the gridless model may be due to the fact 
that it can be derived from the action principle using a truncated Fourier series m and the 
resulting finite dimension dynamical system possesses spectral fidelity to the physical system. The 
advantage of the variational models imiiziiiHi and the spectral fidelity in conservation properties 
and numerical stability deserve further study but will not be pursued here. 

To understand the relative importance of the spectral fidelity in the deposition/interpolation 
and the field solver, in particular the role of amplitude and phase errors, simulations are carried 
out using the momentum conserving PIC model with the deposition/interpolation and the field 
solver implemented in gridless mode respectively. The results are shown in Fig. Operating the 
field solver in the gridless mode effectively turns the finite difference PIC model into a spectral PIC 
model, which improves upon the systematic amplitude error and random spectral errors in a finite 
difference solver. The benefit is seen in Fig. |^where the dominate modes are those separated modes 
(mode 3, 6, 9, 12, 15) up to ujpt — 200, but the onset of whole spectrum excitation, similar to those 
shown in Fig. is delayed to uopt ^ 400 in a longer simulation (not shown). The growth rate 
of mode 15 is about the same as in the case with a finite difference solver and the electron phase 
space shows similar heating. Thus a spectral solver that improves amplitude and random errors 
does not qualitatively change the characteristics of the simulation. Upon adopting the gridless 
deposition and field interpolation (but keeping the finite difference field solver), which improves 
both the systematic phase and amplitude errors, the results are essentially similar to those of the 
gridless model, with a higher floor of the background modes. This confirms that the major loss of 
the spectral fidelity in PIC occurs in the deposition and the field interpolation where their phase 
errors play important roles in FGI. 

The dynamics shown in Fig. [^involves nonlinear coupling of a large range of modes as the system 
evolves, but at the early time the unstable dynamics grows predominately from the interaction 
between density mode k = ±15 and the collective mode iF = 18 which is excited by the initial 
perturbation at iF = 9. To quantitatively determine the behavior of such interaction, we can set 
6x{K^t) = where uo = uor — iuoi, fe(iF,0) = and A = Inserting 

this into the coupled Eqs. (3.11), (3.12) and (3.13), and ignoring the particle aliasing effect then 
A ^ 1, one obtains. 


(kg - ukK) 6 {K - V - i^k'K) 

q, q iVJ 

. ^i(uK+u^,y(u:rt+eg+./2)j^^ {-2kgAoe^^*/Np) {2kq>/Np) , ( 4 . 1 ) 

This can be simplified by equating the time varying phase on both sides of Eq. ( |4.1| ) and setting 
the arguments of the delta functions to zero, i.e., using the following relationships for the mode 
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Figure 3: Simulation mode spectrum from various models. Left column: full mode spectrum ln(|F^(/c)|) as a function 
of time. Right column: evolution of the mode amplitude for the perturbation — mode /c = 9, alias of its second 
harmonic — mode k = 15(= 33 — 9 x 2) and the beat between them — mode k = 6(= 15 — 9). The four rows from 
top to bottom correspond to (1) the momentum conserving (M.C.) PIC model with linear particle shape; (2) the 
momentum conserving PIC model with quadratic particle shape; (3) the energy conserving (E.C.) PIC model with 
quadratic particle shape; and (4) the gridless model with quadratic particle shape. 
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Figure 4: Electron phase space at the end of the simulation (left column), energy (middle column) and momentum 
(right column) conservation from various models. The energy and momentum conservation are shown for ujpAt = 0.2 
(green curves) and cUpAt = 0.02 (blue curves). The three rows from top to bottom correspond to (1) the momentum 
conserving (M.C.) PIC model with quadratic particle shape; (2) the energy conserving (E.C.) PIC model with 
quadratic particle shape; and (3) the gridless model with quadratic particle shape. 
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Figure 5: The mode spectrum (left column), amplitudes of mode 9, 15, 6 (middle column), and electron phase space 
from a PIC model (right column) with (1) a spectral solver (top row); (2) gridless deposition and field interpolation 
(bottom row). 


coupling, 


I'K + I'K' = 1 , 

q = q', 

I'K — kq/K, 

t'K,q^'Z‘, (4.2) 


to give, 


= -2K ^ ^ An^^) J{i-uk) {j'Ke‘^^*An%) . (4.3) 

VK=k^lK J 

where the normalization n n/up, co —>■ Lo/cop, t cOpt, Q —>■ Ng/Np, Q/m —>■ 1, is used and 
An% = 2KAolNp is the normalized initial density perturbation due to the collective mode K. The 
summation in Eq. ( |4.3[ ) is over all integers uk and q that satisfy uk = kq/K = (A: + qkg)/K. For 
the above example, k = ±15, K = 18 and kg = 33, we can rewrite the infinite summation to be. 


E =E 

UK=kqlK leZ 


with g = 6/ =F 1, = 11/ T 1. For 

increases, so in practice only the I = 
coupling is nonlinear. 

Thus, Eq. (|4.3|) becomes. 


\x\ < 0.7, \JiyK{~k'Kx)Ji-UK{k'Kx)\ drops quickly as \uk\ 
0 term needs to be included in the above sum until the 
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(4.4) 

The first (second) term on the right hand side results from the coupling through mode /c = 15 
{k = —15). The coupling from the second term is much stronger than that from the first term 
through mode when the perturbation amplitude is small, i.e., when <C 1. In fact the 

second term leads to absolute instability (cUr = 0), while the first term leads to oscillatory solution 
{(jJi = 0). Therefore, an initially small perturbation grows exponentially until saturation and then 
oscillates at that level. The linear growth rate can be obtained by dropping the first term and 
Taylor expanding the Bessel functions in the second term, which gives. 


-2K 






[-K] 


An%) J 2 + 


s\K 

~\K] 


uji^~s{K)^/\Kj\K]V 

which does not depend on the drift velocity Vq. When 
therefore the normalized instability saturation time is 

taat ~ ln(1.84/An?^)/6c)i. 


(4.5) 


1.84, the RHS of (4.4) is zero, 

(4.6) 

Note that the linear FGI in section 13.21 is different than the case discussed here. For the former 
case, unstable modes grow independently from particle noise. Each unstable mode involves the 
fundamental of the collective motion that can be resolved by the grid (harmonics have been dropped 
in the linear analysis), and the corresponding density perturbation which also contains alias modes 
that cause the instability. For the latter case, the collective motion cannot be resolved by the grid. 
Therefore the density mode contains alias modes only and the instability can be more severe than 
the former case, e.g., for the simulation with large amplitude perturbation. Also, unlike the former 
FGI case which is convective and has stability domains, the latter instability is absolute and has 
no stability domain as long as iF > and s{K) > 0. 

The linear growth rate from Eq. (4.5), uji ~ 0.21, and the instability saturation time from Eq. 


(4.6), tsat ^ 46.5, are compared to the simulation shown in Fig. in which only collective mode 


iF = 18 and electric field modes = ±15 are kept and an initial perturbation in the mode iF = 18 
is applied. Solving Eq. (3.4) for modes k = ±15 gives cji ^ 0.04 which is much smaller than the 


actual grow rate. This simulation result is also in close agreement with the example shown in Fig. 

for the dynamics of mode /c = 15 in the linear growth stage up to the saturation of the instability, 
indicating the early unstable dynamics is due to the coupling between iF = 18 and k = ±15 modes. 
However, the detail dynamics after saturation in the latter simulation involves more mode coupling, 
thus requiring a more complete treatment similar to the above analysis. 


5. Summary 

In this paper, the origin of the FGI is studied by employing particle and spectral resolutions 
into the dynamics of the ID electrostatic PIC model and by contrasting the spectral fidelity of the 
PIC model with respect to the underlying physical system (or the gridless model). The particle 
resolution can be either adopted for individual particle or for the collective motion. The use of 
the particle and spectral resolutions are not the only options for this purpose, but are suitable 
ones as the particle dynamics consist of pair-wise interactions and can also be viewed as collective 
wave-particle interactions. 
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Figure 6: Comparison of the FGI growth (red, dashed, uji ~ 0.21) and the saturation time (black, dotted, tsat ~ 46.5) 
predicted by Eq (4.5) and Eq. ( |4.6| , respectively, with ln\E(k — 15)| in a PIC simulation using a 2nd order particle 
shape and a spectral field solver (green curve). Eq. ( |3.4| gives uoi ^ 0.04. The simulation is initialized with 
a perturbation in the collective mode K — such that An^ = 10“^. Only the collective mode K — IS and the 
electric field modes /c = ±15 are kept in the simulation. Other numerical parameters are the same as in the simulation 
in Eig. [^that employs a spectral solver, no mode filtering and an initial perturbation in mode K = 9 such that 
An^ = 10“^ (blue curve). 


At the individual particle level, the charge deposition and field interpolation operations of the 
PIC models exhibit systematic spectral errors relative to the physical system (or gridless model) 
due to the existence of the spatial alias modes from the use of the discrete grid in conjunction with 
Lagrangian particles in continuous space. In principle, these errors can be calculated for arbitrary 
particle shape, and surprisingly, they have relatively compact analytic forms shown in Eq. (3.1) for 
the B-spline particle shapes and in Eq. (3.2) for the Gaussian particle shape. These forms allow us 
to understand the effects in the spectral domain from using shaped particles in the PIC models. It 
is observed that these systematic spectral errors depend on both kAx and the particle’s normalized 
position in the cell. These errors are directly related to the instability growth rate but they improve 
slowly with the increase of the order (smoothness) of the B-spline particle shape. Another benefit 
of using higher order particle shape is the stronger damping to the short-wavelength modes in the 
fundamental Brillouin zone which can enlarge the stability domain. Such a damping effect also 
applies to the physical modes and may be achieved using a spatial mode filter on the grid which 
has a lower computation cost and does not depend on individual particle’s position. 

At the collective motion level, as it is shown in Eqs. (3.12) and (3.13), the charge deposition and 
field interpolation introduce unphysical mode couplings which eventually result in amplitude and 
phase errors in the collective modes. This is also due to the aliasing of the spatial modes. Within 
the framework of collective wave-particle interactions, Eqs. (3.9), (3.10) and (3.11) are the coupled 
equations describing the evolution of the collective modes in the physical system, while Eqs. (3.12) 
and (3.13) and the discrete time version of Eq. (3.11) are the counterparts for the PIC model, 
whose convergence with respect to time step can be verified to rule out the time discretization as a 
possible origin of the FGI. In fact, early works on FGI already indicated that the numerical heating 
rate is independent of the time step with the leap-frog advance, while our result also shows that the 
growth rates of the unstable modes are not affected by the time step, making it clear the cause of 
FGI can only be the lack of spectral fidelity of Eqs. (3.12) and (3.13) as compared to Eqs. (3.9) and 
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( |3.10 ). In this regard, the difference between the PIC model and the physical system (or gridless 
model) is technical, not fundamental — with a properly-chosen particle shape such that s{k) = 0 
for |A:| > kg/2, we recover Eqs. (3.9) and (3.10) from Eqs. (3.12) and (3.13) for \k\ < kg/2 , thus 
making the PIC models effectively a gridless model. This opens the question of the optimal compact 
particle shapes for the PIC model. Finally, as with many physical instabilities which result from 
the constructive feedback of modes, it can be understood that the systematic spectral errors play 
important roles in the development of FGI. This finding is verified in the simulation comparison 
with the gridless model and may help the design of a PIC model with better stability. 
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Appendix A. Particle shape function and interpolation function 

The function S{r) is used to define particle shape in phase space (for PIC, the shape in velocity 
dimension is assumed to be a Delta function) and the function W{rg, r) is used in the charge 
deposition and field interpolation. The PIC method can be derived from the Vlasov equation 
using the distribution function of shaped particles [ 22 ], which reveals that any symmetric function, 
S{r) — S{—r), or S{k) — S{—k), can be used as long as the averaged fields on the particle is 
defined as in Eq. ( 6 ) of [ 22 | . Self-force consideration and momentum conservation require that the 
same function be used for charge deposition and field interpolation. The other requirement on the 
shape function is J S{r)dr — 1 , or S{k = 0) = 1, while for a valid interpolation function in PIC, 
charge conservation requires r) = 1. From the Poisson summation , this requires 

Yw{rg, r) = = 1, 

rg q 

which is equivalent to IE(fc = 0) = 1 and IE(q • kg) = 0 for n? -h -h ^ 0. Therefore a particle 
shape function is not necessarily an interpolation function, e.g., a linear particle shape function with 
a width of half the grid size is not a valid interpolation function. But valid interpolation functions 
are true particle shape functions, e.g., the family of B-spline interpolation functions used in PIC. 
Furthermore, the convolution of an arbitrary shape function S{r) with these B-spline interpolation 
functions results in valid interpolation functions. In this paper, we will not distinguish the particle 
shape function and the interpolation function unless necessary and will use S{r) for clarity. 

Appendix B. Aliasing and phase factor 

oo _ oo oo 

From the Poisson summation formula ^ ^ikonAx _ ^ S{q — A:oAx/27r) — ^ ~ 

n=—oo q= — oo q= — oo 
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Y^^Uk'-kyvf ^ - kq). 

Vf q 


(B.l) 


We have also used AxAyAz = 1. The infinite sum in Eq. (2.3) is related to the infinite sum on 
the reference grid, Eq. (B.l), by a phase factor. 


^ ^i{k'-k)-Ar, . ^ 87r^^V’(q ' ^rp)6{k' - kq) (B.2) 

rp rf q 

where = '^/ + Avp and '0(q • kg^ ^^p) = is the phase factor. For zero offset 

Avp = (0, 0, 0), '0 = 1; while for a half cell offset, 0 is either 1 or — 1 , e.g. if Avp = (Ax/ 2 , 0, 0), 
0 = (- 1 )^. 


Appendix C. Linear dispersion of electrostatic PIC models 

The linear numerical dispersion of the momentum conserving and energy conserving PIC models 
are given in [1] for finite time step. Here we gives a similar derivation but without introducing 
temporal aliasing, as all PIC quantities, including the distribution function, are discrete in time. 
The perturbed discrete (in time) distribution function is derived in a way similar to Eq. (30) of j8] 
but with one full position update instead of two half updates, 

fi{k,v,cj) = - - -csc[(6l; - k • v)At/2]F{k,cj) • (C.l) 

for which a phase factor is dropped because fi is now needed at the same time step as the 

electric field. 

When the leap-frog update is viewed as a continuous time integration of the Newton’s laws but 
with a force sampled on the time step, temporal aliasing is introduced. In such case the distribution 
function of the eigen mode is defined on the continuous time, and the discrete distribution function 
in PIC is its time sample. The linearization is performed on the former and before the time sampling. 
Using J^ uj-k-v^TT /At ~ ^cot[{uj — k • v)At/2] to evaluate the infinite sum from aliasing [1], one 

arrives at a similar expression to Eq. ( |C.I| ) for the perturbed discrete distribution function but with 
g-z(a;-fc-v)At/ 2 ^g^|-j^^ _ ^ . v)At/2] = cot[( 6 c; — k • v)At/2] — i replaced by cot[( 6 c; — k • v)At/2]. This 
difference does not manifest in the zeroth moment of fi but may appear in higher order moments. 

The normalized equations for the perturbed density pi{k^(jj)^ the potential 0i(fc,6c;),the electric 
field Ei{k^(jj) and the force F{k^v^uj) are. 


pi{k,u) 

— ^^*Sp(fcqr) / fi[kq^V^Uj)dv^ 

q ^ 

(C.2) 

4>iik,Lo) 

= pi{k,uj)/[k]^, 

(C.3) 

Ei{k,u) 

= —iAv(fc) 0 i(fc,Cc;), 

(C.4) 

F{k,u) 

= SE{k)Ei{k,u), 

(C.5) 


All grid quantities are 


where [k]^ is the Poisson operator and K{k) is the gradient operator, 
transformed according to the transform defined in Eq. (2.1). Eq. (C.2) also implies that Avp — 0 
is chosen in Eq. (2.8), so 0(q • kg^ ^^p) = 1- 
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Substituting Eq. (C.l) into the above equations, the dispersion equation is, 


[ \ q 




csc[( 6 c; — kq • v)At/2][K\ • Vvfodv = 0 . (C. 6 ) 


For a cold uniform beam with velocity vq^ it reduces to. 


1 - ^ 


?)■ 


1 

W 


E 


sin^[(6c; — kq • 'Uo)At/2] 


= 0 . 


(C.7) 


Eq. (C.7) has many resonances that are folded into the fundamental Brillouin zone —tt < 
ujAt/2 < TT. For two interacting resonances, the dispersion equation can be simplified as 1 — 
a[sin“^(^) + fein“^(^ + c)] = 0, \b\ < 1. Here we use the ID case and include the funda¬ 
mental mode and the qth {\q\ > 0) alias mode as an example, y = {oj — k • uo)At/2, a = 
At‘^Sp{k)sE{k)Hi{k)k/4:[k]‘^ > 0,b = Sp{kq)sE{kq)E{kq)kq/sp{k)sE{k)E{k)k, and c = {qkgVoAt/2) mod 27r, 
—TV < c < TT. When At is sufficiently small, this equation can be approximated by a quar- 
tic equation in —ac^ -\- 2acy -\- (—a — ab-\- c^) y‘^ — 2cy^ -h = 0 , for which a pair of com¬ 
plex conjugate roots (the other two roots are always real) exist when the discriminant D = 

— 16a^6c^ [(<^(1 + b) — ~ 27a^6c^] < 0. When D > 0, there are (1) four real roots, if a > 0; 

or (2) two pairs of complex conjugate roots, if a < 0. Therefore, for the interaction of any two 
nearby resonances, there are domains of stability/instability defined by the parameters a, 6 and c, 
depending on the signs of D and a. 

Two types of situation can be distinguished. When 6 > 0, the simplified dispersion equation 
resembles that of a physical two-stream instability. Stability requires (a(l + b) — — 27o?bc^ < 0 

and a > 0 , which gives a stable domain 0 < a < ai and unstable domains a > ai or a < 0. Here ai 
is the only real root of (a(l + b) — — 27o?bc^ — 0. For \b\ <C 1, ai (l — 36^/^) c^. When 6 < 0, 

which corresponds to unusual two streams with densities of opposite signs, the stable domain is 
a > ai and the unstable domain is a < ai. 


The dispersion equation Eq. (C.7) applies to both the energy conserving and momentum con¬ 


serving PIC models. The major difference in Eq. (C.7) for these two PIC models is the operator 
K,{kq). For the energy conserving model, K,{kq) = kq^ so b ^ k^ > 0. Therefore, interaction of 
any pair of resonances is of the usual two-stream type and there is a stable domain 0 < a < ai. 
This can be written explicitly, e.g., for the fundamental mode and the qth {\q\ > 0) alias mode, 
[{qkgVoAt/2) mod 2Tr]‘^ > ^Sp{k)sE{k)K{k)k [l - 3{sp{kg)sE{kg)K{kq)kg/sp{k)sE{k)k‘^y/^] t For 
momentum conserving PIC model, E{kq) — kqSmc{kqAx), so b ^ kqSmc{kqAx) which switches sign 
when q switches sign. Therefore, interaction of any pair of resonances can be of either the first or 
the second type. 
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